#ifndef EXTRA_PREC_H
#define EXTRA_PREC_H

#include "/usr/include/qd.h"
#include "/usr/include/x86.h"

#define BEGIN_XPREC unsigned short x86_old_control; x86_fix_start(&x86_old_control);
#define END_XPREC x86_fix_end(&x86_old_control);

typedef dd_real Quad;
typedef qd_real Oct;

const Quad QPI = dd_real::_pi;
const Quad QTWOPI = 2.*QPI;
const Quad QE = dd_real::_e;
const Oct OPI = qd_real::_pi;
const Oct OTWOPI = 2.*OPI;
const Oct OE = qd_real::_e;

inline const Quad& real(const Quad& x) {return x;}
inline const Quad imag(const Quad& x) {return 0.;}
inline const Quad norm(const Quad& x) {return x*x;}
inline const Oct& real(const Oct& x) {return x;}
inline const Oct imag(const Oct& x) {return 0.;}
inline const Oct norm(const Oct& x) {return x*x;}

#include <complex>

inline complex<Quad> operator*(double x,const complex<Quad>& y)
{ return complex<Quad>(x*y.real(),x*y.imag()); }
inline complex<Quad> operator*(const complex<Quad>& x, double y)
{ return y*x; }
inline complex<Quad> operator*=(complex<Quad>& x, double y)
{ x = x*y; return x;}
inline complex<Quad> operator/(double x,const complex<Quad>& y)
{ return x*conj(y)/norm(y); }
inline complex<Quad> operator/(const complex<Quad>& x, double y)
{ return complex<Quad>(x.real()/y,x.imag()/y); }
inline complex<Quad> operator/=(complex<Quad>& x, double y)
{ x = x/y; return x;}

inline complex<Oct> operator*(double x,const complex<Oct>& y)
{ return complex<Oct>(x*y.real(),x*y.imag()); }
inline complex<Oct> operator*(const complex<Oct>& x, double y)
{ return y*x; }
inline complex<Oct> operator*=(complex<Oct>& x, double y)
{ x = x*y; return x;}
inline complex<Oct> operator/(double x,const complex<Oct>& y)
{ return x*conj(y)/norm(y); }
inline complex<Oct> operator/(const complex<Oct>& x, double y)
{ return complex<Oct>(x.real()/y,x.imag()/y); }
inline complex<Oct> operator/=(complex<Oct>& x, double y)
{ x = x/y; return x;}

inline complex<Quad> polar(const Quad& r, const Quad& theta)
{ Quad c,s; sincos(theta,s,c);
    return complex<Quad>(r*c,r*s);
}
inline complex<Oct> polar(const Oct& r, const Oct& theta)
{ Oct c,s; sincos(theta,s,c);
    return complex<Oct>(r*c,r*s);
}

#ifdef GCCCOMP
inline complex<Quad> operator*(const Quad& x,const complex<Quad>& y)
{ return complex<Quad>(x*y.real(),x*y.imag()); }
inline complex<Quad> operator*(const complex<Quad>& x, const Quad& y)
{ return y*x; }
inline complex<Quad> operator*=(complex<Quad>& x, const Quad& y)
{ x = x*y; return x;}
inline complex<Quad> operator/(const Quad& x,const complex<Quad>& y)
{ return x*conj(y)/norm(y); }
inline complex<Quad> operator/(const complex<Quad>& x, const Quad& y)
{ return complex<Quad>(x.real()/y,x.imag()/y); }
inline complex<Quad> operator/=(complex<Quad>& x, const Quad& y)
{ x = x/y; return x;}
inline complex<Oct> operator*(const Oct& x,const complex<Oct>& y)
{ return complex<Oct>(x*y.real(),x*y.imag()); }
inline complex<Oct> operator*(const complex<Oct>& x, const Oct& y)
{ return y*x; }
inline complex<Oct> operator*=(complex<Oct>& x, const Oct& y)
{ x = x*y; return x;}
inline complex<Oct> operator/(const Oct& x,const complex<Oct>& y)
{ return x*conj(y)/norm(y); }
inline complex<Oct> operator/(const complex<Oct>& x, const Oct& y)
{ return complex<Oct>(x.real()/y,x.imag()/y); }
inline complex<Oct> operator/=(complex<Oct>& x, const Oct& y)
{ x = x/y; return x;}
#endif

#endif
